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Abstract 

Bulk  heat  flux  parameterization  is  an  inereasingly  popular  teehnique  for  foreing  non-eoupled  oeean  models.  If  sea  surfaee 
temperature  (SST)  from  the  model  is  eolder  (warmer)  than  observed,  then  the  net  heat  flux  will  be  higher  (lower)  than  observed;  thus, 
bulk  parameterizations  tend  to  keep  model  SST  elose  to  observational  SST  on  long  time  seales.  However,  bulk  parameterizations 
imply  neither  strong  damping  of  SST  variability  nor  strong  relaxation  to  near-surfaee  (e.g.,  at  10  m)  air  temperature  {Tf}.  This  is 
demonstrated  using  SST  simulations  from  a  0.72°  x  0.72°  eos(lat)  (longitude  x  latitude)  resolution  global  HYbrid  Coordinate  Oeean 
Model  (HYCOM)  that  does  not  inelude  assimilation  of  any  SST  data  or  explieit  relaxation  to  any  SST  elimatology,  but  does  use  bulk 
heat  fluxes.  Results  are  diseussed  when  elimatologieal  wind  and  thermal  atmospherie  foreing  for  HYCOM  are  eonstrueted  from  three 
different  arehived  numerieal  weather  predietion  (NWP)  produets:  (1)  the  European  Centre  for  Medium-Range  Weather  Foreeasts 
(ECMWF)  15-year  Re-Analysis  during  1979-1993  (ERA- 15),  (2)  ECMWF  40-year  Re-Analysis  (ERA-40)  during  1979-2002,  and 
(3)  the  Common  Oeean  Referenee  Experiment  Correeted  Normal  Year  foreing  version  1.0  (CORE-CNY)  based  on  the  National 
Centers  for  Environmental  Predietion  (NCEP)  re-analysis  whieh  spans  1948-2002.  To  investigate  the  implieations  of  the  bulk  heat 
flux  approaeh  as  relaxation  to  SST  and  HYCOM  SST  simulations  are  used  to  demonstrate  that  model  SST  errors  with  respeet  to  the 
National  Oeeanie  Atmospherie  Administration  (NOAA)  SST  elimatology  do  not  look  like  SST  fields  from  NWP  produets.  Sueh  a 

result  is  eonfirmed  for  all  simulations  foreed  with  ERA- 15,  ERA-40  and  CORE-CNY,  separately.  Overall,  global  averages  of  mean 
HYCOM  SST  error  are  0.2  °C  (1 .5  °C),  0.4  °C  (1 .7  °C)  and  0.6  °C  (2.3  °C)  with  respeet  to  NOAA  SST  (NWP  T^)  elimatology  when 
the  model  uses  atmospherie  foreing  from  ERA- 15,  ERA-40  and  CORE-CNY,  respeetively. 

©  2008  Elsevier  B.V.  All  rights  reserved. 

Keywords:  Bulk  heat  fluxes;  Ocean  model  SST;  Exchange  coefficients;  Atmospheric  forcing;  Climate 


*  Corresponding  author. 

E-mail  addresses:  alan.wallcraft@nrlssc.navy.mil  (A.J.  Wallcraft), 
birol.kara@nrlssc.navy.mil  (A.B.  Kara), 
harley.hurlburt@nrlssc.navy.mil  (H.E.  Hurlburt), 
echassignet@coaps.fsu.edu  (E.P.  Chassignet), 
ghalliwell@rsmas.miami.edu  (G.H.  Halliwell). 

URL:  http://www7320.nrlssc.navy.mil  (A.J.  Wallcraft). 


1.  Introduction  and  motivation 

Ocean  general  circulation  model  (OGCM)  simula¬ 
tions  are  generally  performed  using  prescribed  atmo¬ 
spheric  forcing  fields,  namely,  momentum  flux  (e.g., 
wind  stress)  and  scalar  forcing  (e.g.,  net  shortwave  and 
longwave  radiation  at  the  sea  surface,  air  temperature 
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and  air  mixing  ratio  at  10  m  above  the  sea  surface). 
These  forcing  fields  are  typically  obtained  from 
archived  NWP  products.  Examples  of  commonly-used 
NWP  products  include  the  European  Centre  for 
Medium-Range  Weather  Forecasts  (ECMWF)  15 -year 
Re-Analysis  (ERA- 15)  (Gibson  et  al.,  1999),  ECMWF 
40-year  Re-Analysis  (ERA-40)  (Kallberg  et  al.,  2004), 
National  Centers  for  Environmental  Prediction  (NCEP) 
Re- Analysis  (Kanamitsu  et  al.,  2002),  and  the  Fleet 
Numerical  Meteorology  and  Oceanography  Center 
(FNMOC)  Navy  Operational  Global  Atmospheric 
Prediction  System  (NOGAPS)  (Rosmond  et  al.,  2002). 

All  NWPs  mentioned  above  use  relatively  sophisti¬ 
cated  boundary  layer  sub-models  to  calculate  sur¬ 
face  fluxes,  but  they  still  depend  on  SST  which  is 
usually  an  analyzed  (non-prognostic)  field.  The  SST 
used  is  now  accurate  enough  that  errors  in  surface  fluxes 
at  the  NWP  grid-scale  are  dominated  by  errors  in 
atmospheric  fields,  such  as  wind  speed  and  cloudiness. 
These  errors  can  be  large,  and  even  the  best  surface 
heat  flux  fields  from  all  sources  (NWP  products  or 
observation-based  climatologies)  do  not  give  a  closed 
global  heat  budget.  In  addition,  using  climatological 
estimates  of  total  heat  flux  to  force  an  ocean  model 
usually  results  in  unrealistic  model  SST,  presenting  an 
inconsistency  with  the  imposed  surface  flux  (Bamier 
et  al.,  1995). 

Even  perfect  fluxes,  with  a  closed  global  heat  budget, 
cannot  be  used  as  the  only  forcing  for  a  stand-alone 
ocean  model  because  the  model’s  “climatology”  is  not 
perfect,  so  such  fluxes  will  lead  to  SST  drift  (e.g., 
Hughes  and  Weaver,  1996).  The  conventional  approach 
to  correcting  this  drift  is  to  add  relaxation  to  observed 
SST,  either  in  place  of,  or  in  addition  to,  prescribed  total 
heat  fluxes  (e.g.,  Bamier  et  al.,  1995;  Maltmd  et  al., 
1998).  As  stated  by  Killworth  et  al.  (2000),  researchers 
at  this  time  did  not  possess  consistent  surface  fields  with 
which  to  force  their  models.  However,  the  use  of  a  bulk 
parameterization  has  become  an  alternative  approach 
because  it  requires  no  explicit  relaxation  to  observed 
oceanic  quantities  (e.g.,  Wallcraft  et  al.,  2003;  Kara 
et  al.,  2003).  In  addition,  using  bulk  parameterizations 
(e.g.,  Parkinson  and  Washington,  1979),  regional 
coupled  ice-ocean  models,  such  as  those  in  the  Baltic 
Sea,  have  been  integrated  over  decades  (Lehmann  and 
Hinrichsen,  2000). 

Bulk  parameterizations  use  the  difference  Ta-SST, 
where  is  the  air  temperature  and  SST  the  ocean  model 
sea  surface  temperature.  The  air  temperature  and  other 
parameters  (air  density  at  the  air-sea  interface,  mixing 
ratio,  and  wind  speed  at  1 0  m  above  the  sea  surface)  are 
typically  from  NWP  products  and  are  used  to  calculate 


the  exchange  coefficients  for  latent  and  sensible  heat 
fluxes  (e.g.,  Kara  et  al.,  2002). 

The  most  important  criticisms  of  bulk  parameteriza¬ 
tions  are  that  (1)  an  infinite  heat  capacity  exists  between 
ocean  and  atmosphere,  i.e.,  it  does  not  give  a  closed 
system,  (2)  since  the  bulk  formula  can  be  linearized  as  a 
function  of  a  difference  between  equilibrium  tempera¬ 
ture  and  SST  (Haney,  1971;  Han,  1984;  Paiva  and 
Chassignet,  2001),  it  can  smear  the  model  SST  when 
the  model  resolution  is  finer  than  the  forcing,  and 
(3)  The  difference,  i.e.,  Ta-SST,  generally  represent 
the  boundary  layer  physics  of  the  NWP  products,  but 
may  be  a  poor  approximation  to  observations.  Note 
that  the  statement  in  (1)  typically  holds  if  the  atmo¬ 
sphere  temperature  is  kept  constant  over  a  longer  period, 
for  example  if  climatological  values  are  used.  Otherwise 
the  atmosphere  temperature  exhibits  a  time  tendency 
too,  which  should  reflect  the  effect  of  the  heat  flux. 

In  this  paper,  we  examine  the  impact  of  bulk  heat  flux 
parameterization  on  the  climatological  SST  produced  by 
simulations  from  a  particular  global  OGCM,  with  no 
assimilation  of,  or  explicit  relaxation  to,  any  observed 
SST  climatology.  In  particular,  reasons  for  preferring  the 
bulk  parameterization  rather  the  net  heat  flux  plus  an 
explicit  relaxation  to  SST  are  discussed. 

A  few  salient  reasons  for  using  SST  from  an  OGCM 
to  discuss  the  bulk  heat  flux  versus  the  direct  relaxa¬ 
tion  approach  are  (1)  SST  plays  an  important  role  in 
air-sea  interactions  through  the  net  heat  flux  at  the 
sea  surface  (e.g.,  Alexander  and  Scott,  1997);  (2)  it  is 
typically  used  in  the  bulk  heat  flux  formulation  to 
parameterize  the  stability  (Kara  et  al.,  2000);  (3)  it 
plays  a  major  role  in  controlling  long  term  climate 
variations,  such  as  the  North  Atlantic  Oscillation 
(Paeth  et  al.,  2003);  (4)  it  is  one  of  the  best  observed 
variables  of  the  upper  ocean  over  the  global  ocean; 
therefore,  SSTs  simulated  from  an  OGCM  can  be 
easily  validated  over  the  globe  (Schopf  and  Loughe, 
1995).  In  addition,  SST  from  an  atmospherically- 
forced  OGCM  (no  assimilation  of  any  ocean  data, 
including  SST)  is  used  because  one  of  our  goals  is  to 
demonstrate  that,  although  the  bulk  heat  flux  formula¬ 
tion  forcing  an  OGCM  surface  temperature  field 
includes  the  knowledge  of  air  temperature  near  the 
sea  surface,  such  use  does  not  imply  an  improper  form 
of  SST  restoring  within  the  model. 

The  paper  is  organized  as  follows.  Section  2 
discusses  the  traditional  bulk  heat  flux  approach  for 
OGCMs.  Section  3  gives  details  of  the  OGCM  used  in 
this  study.  Section  4  presents  SST  results  from  an 
OGCM  in  relation  to  the  bulk  heat  flux  parameteriza¬ 
tion,  and  Section  5  gives  conclusions  of  this  paper. 
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2.  Bulk  heat  flux  approach 

The  total  heat  flux  (2net)  at  the  sea  surface  can  be 
expressed  as  follows: 

Qnet  =  Qsw  ~  Q\w  +  2/  +  2s  5  ( 1 ) 

where  2sw  is  the  net  shortwave  radiation  at  the  sea  surface, 
2iw  is  the  net  longwave  radiation  at  the  sea  surface,  2/ 
is  the  latent  heat  flux,  and  2s  is  the  sensible  heat  flux. 

All  these  individual  components,  and  their  total,  2net? 
are  typically  available  from  the  archived  real-time  and  re¬ 
analysis  data  sets  or  from  climatologies  (ERA-1 5,  ERA-40 
and  CORE-CNY  are  used  here).  Shortwave,  2sw?  and 
longwave,  2iw5  radiation  can  be  directly  measured  at  the 
surface,  but  there  are  far  too  few  such  observations  to 
constrain  even  regional  fluxes.  Thus,  in  practice  all 
components  are  estimated.  All  components  are  affected 
by  surface  type  i.e.,  land  versus  sea  versus  sea-ice  (Garratt 
et  al.,  1 998),  but  here  we  will  only  consider  the  open  ocean. 
Both  shortwave  and  longwave  radiation  depend  heavily  on 
cloudiness,  and  all  flux  components  except  shortwave 
radiation  depend  on  SST  (e.g..  Gill  1982;  Gleckler  and 
Weare,  1997).  HYCOM  uses  a  penetrating  solar  radiation 
scheme  (Kara  et  al.,  2005a)  that  accounts  for  spatial  and 
temporal  water  turbidity  (Kara  et  al.,  2005b,c). 

While  SST  drift  in  OGCMs  may  be  significantly 
reduced  by  relaxation  to  observed  SST  (e.g.,  Bamier 
et  al.,  1995),  such  an  approach  has  the  important  short¬ 
coming  that  the  correct  result  is  prescribed  in  the  model 
simulation.  Rather  than  using  a  direct  SST  relaxation 
scheme,  an  alternative  approach  is  to  calculate  fluxes 
from  a  bulk  parameterization  based  on  near-surface 
atmospheric  fields  using  only  near-surface  wind  speed 
(Va),  or  also  including  air  mixing  ratio  near-surface 
air  temperature  (T^)  (all  of  which  are  10  m  above  the  sea 
surface),  mixing  ratio  for  sea  water  (^s)  and  model  SST, 
as  explained  in  Kara  et  al.  (2005 d). 

Bulk  parameterizations  are  based  on  statistical  fits  to 
observations.  They  are  usually  inexpensive  to  calculate, 
and  the  best  parameterizations  are  generally  very 
accurate.  Examples  of  typical  bulk  parameterizations 
for  latent  and  sensible  heat  fluxes  can  be  found  in 
DeCosmo  et  al.  (1996),  Kara  et  al.  (2000,  2002),  and 
Fairall  et  al.  (2003).  The  formulation  used  in  HYCOM  is 
briefly  described  as  follows: 

a  =  CsCpP,Va(ra  -  SST)  (2) 

Qi  =  CiLp^\^{q^  -  qs)  (3) 

The  air  density  at  the  air-sea  interfaee  (pa)  in  kg 
is  determined  using  the  ideal  gas  law.  The  exchange 


coefficients  (C/  and  Cg)  used  in  HYCOM  include  the 
effects  of  boundary  layer  stability,  and  are  expressed  as 
polynomial  functions  of  Ta-SST,  Va,  and  q^-qs,  the 
latter  through  relative  humidity  (RH),  where  RH  is  at  the 
air-sea  interface  (see  http://www7320.nrlssc.navy.mil/ 
nasec/nasec.html).  In  the  polynomial  functions,  the 
effects  of  water  vapor  flux  in  calculating  the  exchange 
coefficients  are  taken  into  account  through  RH  effects 
that  are  especially  important  at  low  wind  speeds  (Kara 
et  al.,  2005d).  The  total  heat  flux  (Eq.  (1))  is  therefore 
expressed  as  a  polynomial  function  of  the  model  SST 
and  can  be  linearized  to  first  order  as 

2net  =  ^(2^*-SST),  (4) 

in  which  the  apparent  equilibrium  temperature  and 
the  relaxation  coefficient  A  are  time  and  space 
dependent  and  can  be  computed  from  the  atmospheric 
fields  (Haney,  1971;  Han,  1984;  Paiva  and  Chassignet, 
2001).  We  have  included  2sw  in  the  Haney  parameter¬ 
ization,  as  is  usually  done,  even  though  it  does  not 
depend  on  SST.  The  total  heat  flux  therefore  implies  a 
temperature  relaxation,  but  not  toward  the  atmospheric 
temperature  Ta,  nor  toward  a  “correcf’  SST  (Tc).  The 
equilibrium  temperature  indeed  must  differ  from  T^, 
as  it  would  otherwise  imply  zero  heat  flux  when  the 
model  SST  is  equal  to  T^. 

2.1.  Preference  for  a  bulk  parameterization 

Instead  of  using  a  direct  relaxation  to  any  SST 
climatology,  a  common  alternative  has  been  to  use  NWP 
net  flux  plus  an  explicit  relaxation  to  the  correct  SST, 
Since  SST  is  one  of  the  best  observed  geophysical  fields, 
why  should  a  bulk  parameterization  be  preferred  over 
the  NWP  net  flux  plus  an  explicit  relaxation  to  SST? 
One  answer  is  that  the  e-folding  time  implicit  in  accurate 
bulk  parameterizations  is  representative  of  the  air-sea 
interface.  It  is  possible  to  construct  explicit  relaxation 
terms  that  are  patterned  on  bulk  parameterizations,  or 
alternatively  to  use  measured  climatological  e-folding 
times  from  observations  or  NWP  fields  (or  the  results  of 
a  bulk  parameterization  applied  to  NWP  fields).  These 
approaches  are  all  of  the  Haney  type  (Haney,  1971 ;  Han, 
1984;  da  Silva  et  al.,  1994;  Paiva  and  Chassignet,  2001) 
as  follows 

e(ssT)  =  e(rs)  +  |p|rjre-ssT),  (5) 

where  the  effective  SST  relaxation  e-folding  time  scale 
implied  by  |  depends  on  many  factors.  However,  it 
is  certainly  positive  and  typically  large  enough  that 
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model  SST  does  not  exhibit  long  term  drift.  In  practice 
relaxation  terms  are  often  simpler  than  this,  i.e.,  a 
constant  e-folding  time  and  some  approximation  to  T^. 

Eddy-resolving,  atmospherically-forced  ocean  models 
without  ocean  data  assimilation  should  not  relax  strongly 
to  synoptic  observed  SST  because  the  observed  SST  will 
include  SST  anomalies  associated  with  eddies  and 
western  boundary  currents  at  different  locations  than  in 
the  model  simulation.  This  is  because  of  flow  instabilities 
in  the  model  which  are  not  a  deterministic  response  to  the 
forcing.  The  bulk  parameterization  does  not  use  the  NWP 
SST  (Ts)  explicitly,  and  the  NWP  fields  are  often  on  a 
much  coarser  grid  than  the  ocean  model  one  that  barely 
resolves  such  details  as  eddy  and  ocean  current  locations. 
However,  this  raises  the  question  of  to  what  degree  the 
implied  relaxation  in  the  heat  flux  can  smear  the  model 
SST  when  the  model  resolution  is  finer  than  the  forcing.  It 
also  raises  the  question  as  to  whether  the  NWP  atmo¬ 
spheric  temperature  is  strongly  correlated  to  the  NWP 
SST  Ts,  therefore  introducing  oceanic  information  in  the 
heat  flux  formulation.  In  contrast,  eddy-resolving  stand¬ 
alone  ocean  models  that  assimilate  sea  surface  height 
(from  satellite  altimeters)  have  ocean  fronts  and  eddies  in 
the  observed  location  (Smedstad  et  al.,  2003;  Chassignet 
et  al.,  2006;  Shriver  et  al.,  2007)  and  therefore  can  relax 
directly  to  high  resolution  observed  SST,  although  they 
may  instead  (or  in  addition)  directly  assimilate  SST.  From 
the  preceding  discussion,  simple  relaxation  is  never 
appropriate  as  the  only  heat  flux  term,  so  it  must  be 
augmented  either  with  the  NWP  net  heat  flux  or  a  bulk 
parameterization  heat  flux. 

Regions  where  the  ocean  model’s  mean  SSTs  are  less 
accurate  can  often  be  traced  to  deficiencies  in  the 
atmospheric  fields.  In  particular,  systematic  biases  in 
wind  speed  will  always  lead  to  poor  fluxes  from  the  bulk 
parameterization  and  biases  in  cloudiness  will  have  a  large 
effect  on  shortwave  radiation  (partially  offset  by  changes 
to  downward  longwave  radiation).  Biased  winds  are 
unlikely  to  provide  better  surface  heat  fluxes  from  the 
sophisticated  boundary  layer  sub-model  in  the  NWP  than 
from  the  bulk  parameterization.  Assimilation  of  wind 
speeds  from  satellites  should  improve  the  accuracy  of 
NWP  surface  winds,  but  there  are  still  large  biases  in  some 
coastal  regions,  and  there  may  be  smaller  systematic 
biases  on  larger  scales.  Cloud  cover  is  difficult  to  obtain 
accurately,  and  known  to  be  deficient  regionally  in  both 
archived  real-time  and  re-analysis  products. 


the  SST,  therefore  introducing  oceanic  information  into 
the  heat  flux  formulation,  the  long  term  (climatological) 
mean  of  the  difference  between  near-surface  and  SST 
from  an  observation-based  data  set  and  the  SST  used  by 
NWP  products  are  investigated  over  the  global  ocean 
(Fig.  1).  Specifically,  the  observation-based  data  set  is 
the  Comprehensive  Ocean  Atmosphere  Data  Set 


Climatological  mean  bias:  Ta-Ts 


2.2.  Relationship  between  air  temperature  and  SST 

In  order  to  answer  the  question  as  to  whether  the 
atmospheric  temperature,  Ta,  is  strongly  correlated  to 


Fig.  1.  Climatological  mean  of  the  differenee  between  near-surface  air 
temperature  (T^)  and  sea  surface  temperature  (i.e.,  T^-T^  over  the 
global  ocean.  Both  and  f  fields  from  COADS,  ERA- 15,  ERA-40 
and  CORE-CNY  are  interpolated  to  the  0.72°  HYCOM  grid,  and  the 
difference  between  the  two  is  formed  at  each  grid  point. 
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(COADS)  climatology  based  largely  on  ship  observa¬ 
tions  and  buoy  measurements  during  1945-1989  (da 
Silva  et  al.,  1994)  and  the  NWP  data  products  are  the  re¬ 
analysis  products,  ERA-15  (1979-1993),  ERA-40 
(1979-2002)  and  CORE-CNY  (1948-2002).  COADS 
is  the  only  observation-based  climatology  among  these 
products,  and  is  intended  to  complement  the  compar¬ 
isons.  This  is  the  new  1/2°  x  1/2°  climatology  based  on 
the  atlas  of  surface  marine  data.  Supplement  B  (http:// 
www.nodc.noaa.gov/OC5/asmdnew.html).  Note  that 
hereinafter  the  COADS  SST,  as  well  as  the  SSTs  used 
in  the  NWP  products  will  be  denoted  as  T^. 

Although  the  time  periods  over  which  each  climatol¬ 
ogy  (COADS,  ERA- 15,  ERA-40  and  NCEP)  is 
constructed  are  different,  they  are  all  >15  years  which 
is  adequate  to  represent  a  long  term  mean  over  the 
global  ocean.  Near-surface  and  data  for  each 
product  are  available  from  the  National  Center  Atmo¬ 
spheric  Research  (NCAR)  web  site  (available  online  at 
http://dss.ucar.edu/catalogs/).  COADS  provides 
monthly  mean  fields  directly,  while  for  the 

three  NWP  products  we  form  climatological  means 
using  Ta  and  the  fields  archived  at  sub-daily  intervals. 
Near-surface  atmospheric  temperature  from  ERA- 15, 
ERA-40  and  NCEP  is  at  2  m,  while  that  from  COADS  is 
at  10  m  above  the  sea  surface.  The  adjustment  from  2  to 
10  m  Ta  is  very  small  as  indicated  by  the  comparison  of 
Ta  values  in  Table  1  and  hence  was  ignored. 

When  sea-ice  is  present,  surface  is  over  sea-ice  in 
the  re-analysis  products,  but  not  in  COADS.  The  field 
from  NCEP/NCAR  is  actually  obtained  from  the 
CORE-CNY  (Corrected  Normal  Year)  data  set  which 

Table  1 

Global  averages  for  near-surface  air  temperature  (Ta),  f  and  difference 
between  the  two  from  the  observational-based  COADS  data  set  and 
three  NWP-based  products 


Data 

T 

T 

T  -T 

COADS 

21.4 

22.0 

-0.6 

ERA- 15 

21.4 

22.1 

-0.7 

ERA-40 

21.2 

22.1 

-0.9 

COARE-CNY 

20.9 

22.0 

-1.1 

Further  information  about  COADS,  ERA- 15,  ERA-40  and  CORE- 
CNY  can  be  found  online  at  several  web  sites.  For  example,  as  of  this 
writing,  the  web  addresses  are  http://iridl.ldeo.columbia.edu/ 
SOURCES/.DASILVA/.SMD94  for  COADS,  http://badc.nerc.ac.uk/ 
data/ecmwf-era/ERA.html  for  ERA- 15,  http://badc.nerc.ac.uk/data/ 
ecmwf-e40/e40_background_html  for  ERA-40,  and  http://datal.gfdl. 
noaa.gOv/nomads/forms/mom4/CORE/CNYF_lpO.html  for  CORE- 
CNY.  The  web  addresses  may  be  subject  to  changes  by  the  originators. 
All  the  NWP  products  have  different  boundary  layer  parameteriza- 
tions,  physics,  data  assimilation  methods  and  different  satellite  data 
used  in  the  assimilations.  Therefore,  differences  in  their  output 
variables  are  expected. 


corrected  for  excessively  cold  values  in  the  Antarctic. 
Ignoring  high  latitudes  where  sea-ice  forms,  the  fields 
are  broadly  similar  to  each  other  with  climatological 
warmer  than  nearly  everywhere  but  usually  by 
<  1  °C.  The  exception  is  in  areas  where  strong  currents 
transport  water  which  is  significantly  warmer  than  the 
air  in  the  mean,  e.g..  Gulf  Stream,  Kuroshio  and 
Agulhas.  Further  details  can  be  found  in  Kara  et  al. 
(2007). 

The  same  long  term  mean  fields  shown  in  Fig.  1  are 
plotted  as  scatter  diagrams  of  versus  over  the 
global  ocean  (Fig.  2).  The  regions  where  ice  is  located 
are  excluded  in  the  evaluation  procedure.  There  is  a 
strong  linear  relationship  between  and  SST  with  R 
values  >0.99  for  all  data  sets.  This  is  clearly  evident 
from  the  slope  values  being  close  to  1  in  all  cases  (1.008, 
0.991,  1.012  and  1.033  for  COADS,  ERA-15,  ERA-40 
and  CORE-CNY,  respectively).  Intercept  values  are  also 
similar  with  values  of  0.468,  0.911,  0.711  and  0.405. 
The  difference  between  and  (i.e.,  T^-Tf)  ranges 
from  -  0.64  °C  for  COADS  to  - 1 .09  °C  for  CORE-CNY 
(Table  1). 

In  contrast,  near-surface  and  T^-T^  are  only 
weakly  correlated  (Fig.  3)  with  linear  correlation  coef¬ 
ficients  of-0.12,  -0.14,  -0.17  and -0.25  for  COADS, 
ERA- 15,  ERA-40  and  CORE-CNY.  None  of  the  cor¬ 
relation  values  are  statistically  significant  in  comparison 
to  a  correlation  value  of  0  at  a  95%  confidence  interval 
based  on  the  student’s  t-tQSt.  Thus,  itself  does  not 
have  strong  influence  on  T^-T^. 

3.  HYCOM  description 

The  HYbrid  Coordinate  Ocean  Model  (HYCOM)  is 
based  on  a  primitive-equation  formulation  discussed  by 
Bleck  (2002),  in  detail.  There  have  been  several 
HYCOM  applications  investigating  a  variety  of  pro¬ 
cesses  in  different  ocean  basins  and  enclosed  seas. 
Examples  of  such  studies  include  the  North  Atlantic 
(Chassignet  et  al.,  2003;  Halliwell  2004;  Thacker  et  al., 
2004),  the  Indian  Ocean  (Han  et  al.,  2004),  the  tropical 
Pacific  (Shaji  et  al.,  2005),  and  the  Black  Sea  (Kara 
et  al.,  2005a,b,c).  HYCOM  has  also  been  used  as  the 
ocean  component  of  a  coupled  atmosphere-ocean 
model  in  global  climate  studies  (e.g..  Sun  and  Hansen 
2003;  Bleck  and  Sun  2004),  and  for  eddy-fresolving 
ocean  prediction  (Chassignet  et  al.,  2006). 

For  the  present  study,  we  will  introduce  HYCOM 
configured  for  the  global  ocean,  including  the  latest 
advances  in  the  model  development  (see  Appendix  for 
details).  The  model  presented  here  is  a  stand-alone 
ocean  model  with  no  assimilation  of  any  ocean  data. 
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including  SST,  and  no  relaxation  to  any  other  data 
except  sea  surface  salinity  (SSS)  to  keep  the  evaporation 
minus  precipitation  balance  on  track.  General  features  of 
the  global  atmospherically-forced  HYCOM  are  given  in 
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Fig.  2.  Scatter  plots  of  versus  f  based  on  annual  mean  values,  each 
at  1°  bins  over  the  global  ocean.  Both  and  f  fields  are  from 
COADS,  ERA- 15,  ERA-40  and  CORE-CNY,  respectively. 


Fig.  3.  Scatter  plots  of  near-surface  air  temperature  {T^  versus 
based  on  annual  mean  values  (see  Fig.  1),  each  at  1°  bins  over  the 
global  ocean. 
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Section  3.1  and  the  atmospheric  forcing  used  in  the 
model  simulations  is  explained  in  Section  3.2. 

3.1.  General  features  of  global  HYCOM 

HYCOM  contains  five  prognostic  equations:  two  for 
the  horizontal  velocity  components,  a  mass  continuity  or 
layer  thickness  tendency  equation  and  two  conservation 
equations  for  a  pair  of  thermodynamic  variables,  such  as 
salt  and  potential  temperature  or  salt  and  potential 
density  (Bleck,  2002).  The  model  behaves  like  a 
conventional  cr  (terrain-following)  model  in  very 
shallow  oceanic  regions,  like  a  z-level  (fixed-depth) 
coordinate  model  in  the  mixed  layer  or  other  unstratified 
regions,  and  like  an  isopycnic-coordinate  model  in 
stratified  regions.  However,  the  model  is  not  limited  to 
these  coordinate  types  (Chassignet  et  al.,  2003). 
Typically,  HYCOM  has  isopycnal  coordinates  in  the 
stratified  ocean  but  uses  the  layered  continuity  equation 
to  make  a  dynamically  smooth  transition  to  z-levels  in 
the  unstratified  surface  mixed  layer  and  cr-levels  in 
shallow  water.  The  optimal  coordinate  is  chosen  every 
time  step  using  a  hybrid  coordinate  generator.  Thus,  the 
model  automatically  generates  the  lighter  isopycnal 
layers  needed  for  the  pycnocline  during  summer,  but 
during  winter  the  same  layers  define  fixed-depth  z-level 
coordinates. 

HYCOM  domain  used  in  this  paper  spans  the  global 
ocean  from  78°S  to  90°N.  It  has  a  0.72°  equatorial 
Mercator  grid  between  78°S-47°N,  with  an  Arctic  bi¬ 
polar  patch  above  47°N.  Average  zonal  (longitudinal) 
grid  resolution  for  the  0.72°  global  model  varies  from 
^  80  km  at  the  equator  to  ~  60  km  at  mid-latitudes  (e.g., 
at  40°  N).  The  meridional  (latitudinal)  grid  resolution  is 
doubled  near  the  equator  to  better  resolve  the  equatorial 
wave-quide  and  halved  in  the  Antarctic  for  computa¬ 
tional  efficiency.  Hereinafter,  the  model  resolution  will 
be  referred  to  as  0.72°  for  simplicity.  Zonal  and 
meridional  array  lengths  are  500  and  457,  respectively. 
At  this  resolution,  coastal  regions  are  not  represented  in 
great  detail  (so  sigma-levels  are  not  used),  and  the  model 
land-sea  boundary  is  at  the  50  m  isobath  (with  a  closed 
Bering  Strait). 

All  global  simulations  use  sigmaO  (sigma-theta),  i.e. 
potential  density  referenced  to  the  surface  (0  dbar)  with 
no  thermobaric  correction.  There  are  26  hybrid  layers  in 
the  vertical.  In  general,  the  model  needs  fewer  vertical 
coordinate  surfaces  than,  say,  a  conventional  z-level 
model,  because  isopycnals  are  more  efficient  in 
representing  the  stratified  ocean,  as  discussed  in  Hurlburt 
et  al.  (1996)  and  Kara  et  al.  (2005a).  The  target  density 
values  for  the  isopycnals  and  the  decreasing  change  in 


density  with  depth  between  isopycnal  coordinate 
surfaces  are  based  on  the  1/4°  Generalized  Digital 
Environmental  Model  (GDEM)  climatology  (NAVO- 
CEANO,  2003).  The  density  difference  values  were 
chosen  so  that  the  layers  tend  to  become  thicker  with 
increasing  depth,  with  the  lowest  abyssal  layer  being  the 
thickest.  The  near-surface  z-level  regime  is  a  natural 
consequence  of  HYCOM ’s  minimum  layer  thickness.  In 
this  case,  the  minimum  thickness  of  layer  1  is  3  m,  and 
this  minimum  increases  1.125x  per  layer  up  to  a 
maximum  at  12  m,  and  target  densities  are  chosen, 
such  that  at  least  the  top  four  layers  are  in  z-level 
coordinates. 

The  simulations  use  realistic  bottom  topography 
constructed  from  the  NRL  2  min  resolution  bathymetry. 
Extensive  quality  control  of  bottom  topography  was 
performed  in  straits  and  near  coastlines.  Monthly  mean 
temperature  and  salinity  from  the  GDEM  climatology  in 
August  are  used  to  initialize  the  model.  There  is  a 
relaxation  to  monthly  mean  SSS  from  the  Polar  science 
center  Hydrographic  Climatology  (PHC).  The  PHC 
climatology  is  chosen  for  its  accuracy  in  the  Arctic 
region  (Steele  et  al.,  2001).  The  reference  mixed  layer 
thickness  for  the  SSS  relaxation  is  30  m  (30  days  in  30  m 
e-folding  time).  The  actual  e-folding  time  depends  on 
the  mixed  layer  depth  (MED)  and  is  30  x  30/MLD  days, 
i.e.  it  is  more  rapid  when  the  MED  is  shallow  and  less 
so  when  it  is  deep.  Such  a  relaxation  is  necessary  to 
prevent  SSS  drift,  and  is  in  addition  to  the  evaporation- 
precipitation  budget. 

HYCOM  has  several  mixed  layer/vertical  mixing 
options  (see  Halliwell  (2004)  for  a  discussion  and 
evaluation).  In  this  paper,  the  non-slab  K-Profile 
Parameterization  (KPP)  of  Large  et  al.  (1997)  is  used. 

3.2.  Atmospheric  forcing 

The  model  reads  in  the  following  time-varying 
atmospheric  forcing  fields:  wind  forcing  (zonal  and 
meridional  components  of  wind  stress,  wind  speed  at 
10  m  above  the  sea  surface)  and  thermal  forcing  (air 
temperature  and  air  mixing  ratio  at  10  m  above  the  sea 
surface,  precipitation,  net  shortwave  radiation  and  net 
longwave  radiation  at  the  sea  surface).  The  wind/ 
thermal  forcing  was  constructed  from  three  NWP 
products:  (i)  1.125°  x  1.125°  ERA-15,  (ii)  1.125°  x 
1.125°  ERA-40  and  (hi)  1.875°  x  1.875°  CORE-CNY. 
When  using  each  product,  the  goal  is  to  include  high 
temporal  variability  in  the  forcing,  while  maintaining  a 
climatology. 

The  original  atmospheric  forcing  data  set  from  ERA-1 5 
spans  the  period  1  January  1979-31  December  1993. 
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The  climatological  ERA- 15  forcing  used  here  consists  of 
15 -year  monthly  averages  over  1979-1993,  interpolated  in 
time  to  6  hourly  and  with  6  hourly  sub-monthly  anomalies 
from  operational  ECMWF  in  September  1994  to  Septem¬ 
ber  1995  added  to  the  winds  only  Choosing  another  time 
period  for  the  6  hourly  wind  anomalies  (other  than  1994- 
95)  did  not  make  any  significant  impact  on  the  model  SST. 
ECMWF  used  a  spectral  model  to  generate  the  ERA- 15 
dataset.  The  ERA- 15  re-analysis  project  incorporated  a 
number  of  in-situ  and  satellite-based  data  types  over  the  1 5- 
year  period.  SST  analyses  for  ERA- 15  are  provided  by  the 
United  Kingdom  Meteorological  Office  (UKMO)  for  the 
early  period,  and  by  NOAA  from  November  1981 
onwards.  Sea-ice  cover  has  been  derived  from  satellite  data. 

The  ERA-40  project  applies  a  modem  variational 
data  assimilation  technique  for  the  past  conventional 
and  satellite  observations.  The  model  physics  and  the 
surface  parameterization  have  been  upgraded  and 
improved  since  ERA- 15.  A  significant  difference 
between  ERA-40  and  ERA- 15  is  in  the  use  of  satellite 
data.  ERA-40  uses  the  Advanced  Tiros  Operational 
Vertical  Sounde  (ATOVS)  radiances  directly,  while  in 
ERA- 15  temperature  and  humidity  retrievals  were  used. 
In  addition.  Special  Sensor  Microwave  Imager  (SSM/I) 
and  European  Remote  Sensing  Satellite  (ERS)  data  are 
used  in  ERA-40.  The  SST/Ice  data  set  produced  by  the 
Hadley  Centre  and  National  Oceanic  and  Atmospheric 
Administration  (NOAA)/National  Environmental  Satel¬ 
lite,  Data,  and  Information  Service  (NESDIS)  has  been 
made  available  to  the  ERA-40  project.  For  use  in 
the  analyses  and  HYCOM  simulations  performed  here, 
25-year  monthly  climatologies  of  atmospheric  forcing 
parameters  were  formed  from  the  ERA-40  archives 
spanning  1979-2002  with  the  same  6  hourly  wind 
anomalies  added  as  in  the  ERA- 15  case. 

The  atmospheric  forcing  from  CORE-CNY  is  for  a 
single  year  and  combines  re-analysis  data  from  NCEP 
with  satellite  measurements  to  reduce  errors  existing  in 
the  original  NCEP  fields,  especially  radiation  fields 
(e.g.,  Lee  et  al.,  2005).  Heat  fluxes  for  the  CORE-CNY 
forcing  are  produced  by  applying  the  NCAR  bulk 
parameterization  to  NOAA  SSTs,  however  HYCOM 
calculates  its  own  latent  and  sensible  heat  fluxes  based 
on  model  SST  (see  below).  The  forcing  includes  6 
hourly  representative  variability  in  all  its  fields. 

In  addition  to  wind  and  thermal  forcing,  HYCOM 
forcing  incorporates  monthly  mean  climatologies  of  river 
discharge  and  satellite-based  attenuation  coefficients  for 
Photosynthetically  Active  Radiation  (A:par  in  m~^).  The 
shortwave  radiation  at  depth  is  calculated  using  a  spatially 
and  temporally  varying  monthly  A:par  climatology  as 
processed  from  the  daily-averaged  (attenuation 


coefficient  at  490  nm)  data  set  from  Sea-viewing  Wide 
Field-of-view  Sensor  (SeaWiFS)  during  1997-2001. 
Thus,  using  ocean  color  data,  the  effects  of  water  turbidity 
are  included  in  the  model  simulations  through  the 
attenuation  depth  in  m)  for  shortwave  radiation. 

The  rate  of  heating/cooling  of  model  layers  in  the  upper 
ocean  is  obtained  from  the  net  heat  flux  absorbed  from  the 
sea  surface  down  to  a  depth,  including  water  turbidity 
effects  (e.g.,  Kara  et  al.,  2005a).  Previously,  it  was  shown 
that  water  turbidity  can  be  quite  significant  in  SST 
simulations,  especially  in  equatorial  regions  (e.g.,  Kara 
et  al.,  2004,  2005b). 

The  net  surface  heat  flux  that  has  been  absorbed  (or 
lost)  by  the  upper  ocean  to  depth  is  parameterized  as  the 
sum  of  the  downward  surface  solar  radiance,  upward 
longwave  radiation,  and  the  downward  latent  and 
sensible  heat  fluxes  (see  Section  2).  Net  solar  radiation 
(the  sum  of  net  shortwave  and  longwave  radiation)  at  the 
sea  surface  is  so  dependent  on  cloudiness  that  it  is  taken 
directly  from  the  given  NWP  product  (ECMWF  or 
CORE-CNY)  for  use  in  the  HYCOM.  The  net  longwave 
flux  is  the  sum  of  downward  longwave  (from  the 
atmosphere)  and  upward  black-body  radiation.  The 
NWP  (input)  black-body  radiation  is  corrected  within 
HYCOM  to  allow  for  the  difference  between  NWP  SST, 
T;,  and  HYCOM  SST  (Kara  et  al.,  2005b).  The 
downward  longwave  radiation  is  often  not  archived, 
but  if  archived,  input  net  longwave  is  calculated  using 
an  archived  NWP  SST. 

Latent  and  sensible  heat  fluxes  at  the  air-sea  interface 
are  calculated  using  computationally  efficient  bulk 
formulae  that  include  the  effects  of  dynamic  stability 
(Kara  et  al.,  2005d).  The  details  of  the  parameterizations 
are  given  in  Section  2.  HYCOM  treats  rivers  as  a  “runoff’ 
addition  to  the  surface  precipitation  field.  The  flow  is  first 
applied  to  a  single  ocean  grid  point  and  smoothed  over 
surrounding  ocean  grid  points,  yielding  a  contribution  to 
precipitation  in  m  s~  ^ 

4.  HYCOM  SST  simulations 

In  this  section,  we  investigate  whether  or  not  monthly 
mean  SSTs  obtained  from  atmospherically-forced 
HYCOM,  which  includes  the  effects  of  bulk  parameter¬ 
izations  in  its  surface  energy  balance  (see  Section  2),  is 
strongly  controlled  by  the  near-surface  used  in  the 
atmospheric  forcing.  The  global  HYCOM  simulations 
presented  in  this  paper  were  performed  with  no  assimila¬ 
tion  of  any  oceanic  data  except  initialization  from 
climatology.  There  is  only  weak  relaxation  to  sea  surface 
salinity  to  keep  the  evaporation-precipitation  budget  on 
track  in  the  model.  Model  simulations  are  performed 
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using  atmospheric  wind/thermal  forcing  from  each  NWP 
product  (ERA- 15,  ERA-40  and  CORE-CNY),  separately. 
We  did  not  perform  any  simulations  using  the  alternative 
approach  of  applying  NWP  net  heat  fluxes  and  relaxing  to 
“observed”  SST,  both  because  this  explicitly  includes  the 
“answer”,  at  least  for  SST,  in  the  model  run  and  because 
we  would  have  to  choose  one  particular  relaxation  scheme 
which  would  still  leave  open  the  question  of  whether 
some  other  variant  would  be  more  appropriate. 

Near  statistical  equilibrium  was  reached  for  each 
simulation  after  four  model  years.  A  linear  regression 
analysis  was  performed  for  domain  averaged  quantities 
(layer  temperature,  salinity,  potential  and  kinetic  energy, 
etc.)  to  investigate  statistical  equilibrium  in  each  layer, 
and  is  expressed  numerically  as  %  change  per  decade. 
The  model  is  deemed  to  be  in  statistical  equilibrium 
when  the  rate  of  potential  energy  change  is  acceptably 
small  (e.g.,  <  1%  in  5  years)  in  all  layers.  The  statistical 
equilibrium  was  accomplished  using  climatological 
monthly  mean  thermal  atmospheric  forcing,  but  with 
wind  forcing  that  includes  the  6-h  sub-monthly  varia¬ 
bility  because  of  mixed  layer  sensitivity  to  high 
frequency  forcing  (e.g.,  Kara  et  al.,  2005a). 

Performing  a  1-year  global  HYCOM  simulation  using 
any  of  the  atmospheric  forcing  products  required  ~  1 5 
wall-clock  hours  on  64  HP/Compaq  SC45  processors. 
Thus,  the  0.72°  global  HYCOM  provides  inexpensive, 
non-eddy-resolving  ocean  model  simulations. 

Monthly  mean  SST  fields  obtained  from  the  model 
simulations  are  evaluated  through  extensive  model-data 
comparisons  using  various  statistical  metrics  (Section 
4.1).  Within  the  quantitative  framework,  statistical  error 
and  skill  analyses  are  then  presented  for  quality 
assessment  of  the  model  in  simulating  climatological 
monthly  mean  SST  using  bulk  parameterizations  and  the 
three  atmospheric  forcing  products  (Section  4.2). 


are  mainly  designed  for  large-scale  climate  studies. 
They  are  derived  by  a  linear  interpolation  of  the  weekly 
optimal  interpolation  (01),  and  use  in-situ  and  satellite 
SSTs,  making  it  a  reliable  candidate  for  model-data 
comparisons  over  the  global  ocean.  The  existence  of  the 
ice  field  in  the  NO  A  A  data  set  is  also  an  advantage  for 
the  model  SST  validation  in  the  Arctic  and  Antarctic. 
We  will  also  use  climatologies  from  NWP  products  to 
compare  with  HYCOM  SSTs  in  Section  4.2. 

Monthly  mean  HYCOM  SST  obtained  from 
HYCOM  simulations  using  wind  and  thermal  forcing 
from  ERA- 15,  ERA-40  and  CORE-CNY  is  validated 
using  mean  error  (ME),  root-mean-square  (RMS)  SST 
difference,  correlation  coefficient  {R)  and  non-dimen¬ 
sional  skill  score  (SS).  Let  Y,  (/=  1,  2,  -,  12)  be  the  set  of 
monthly  mean  NOAA  reference  (observed)  SST  values 
from  January  to  December,  and  let  7^  (i=  1,  2,  -,  12)  be 
the  set  of  corresponding  HYCOM  estimates  at  a  model 
grid  point.  Also  let  Y  (7)  and  Gx  (cry)  be  the  mean  and 
standard  deviations  of  the  reference  (estimate)  values, 
respectively.  The  statistical  relationships  (e.g.,  Murphy, 
1995)  between  NOAA  and  HYCOM  SST  time  series  at 
a  given  grid  point  are  expressed  as  follows: 


ME  =  7-Y 


RMS  = 
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4.1.  Statistical  metrics 

Monthly  mean  HYCOM  SST  climatologies  are 
constructed  from  SSTs  obtained  from  the  last  4  (out  of 
8)  model  years.  For  example,  the  climatological  mean 
January  SST  is  formed  using  mean  January  SSTs  from 
model  years  5-8.  The  same  process  is  repeated  for  the 
other  months.  Because  all  simulations  performed  with 
the  0.72°  model  use  climatological  mean  forcing  from 
NWP  products  and  there  are  no  mid-latitude  mesoscale 
eddies,  the  4-year  averaging  period  is  sufficient. 

For  model-data  comparisons  the  NOAA  SST  clima¬ 
tology  (Reynolds  et  al.,  2002)  is  taken  as  a  reference 
(truth)  because  its  resolution  (1°  x  1°)  is  close  to  that  of 
HYCOM  (0.72°  X  0.72°  cos(lat)).  The  NOAA  SST  fields 


ME  (i.e.,  annual  bias)  is  the  mean  error  between  the 
mean  HYCOM  and  NOAA  SST  values,  RMS  (root- 
mean-square)  SST  difference  is  an  absolute  measure  of 
the  distance  between  the  two  time  series,  and  the  R 
value  is  a  measure  of  the  degree  of  linear  association 
between  the  time  series.  The  non-dimensional  SS  given 
in  Eq.  (9)  is  the  fraction  of  variance  explained  by 
HYCOM  minus  two  dimensionless  biases  (conditional 
bias,  ^cond?  and  unconditional  bias,  ^uncond)  which  are 
not  taken  into  account  in  the  R  formulation  (8)  as 
explained  in  Murphy  (1988),  in  detail,  ^uncond?  described 
as  systematic  bias,  is  a  measure  of  the  difference  between 
the  means  of  NOAA  and  HYCOM  SST  time  series,  ^cond 
is  a  measure  of  the  relative  amplitude  of  the  variability  in 
the  NOAA  and  HYCOM  time  series  or  simply  a  bias  due 
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to  differences  in  standard  deviations  of  the  SST  time 
series.  Note  SS  is  1.0  (negative)  for  perfect  (poor) 
HYCOM  SST  simulations. 

4.2.  HYCOM  SST  evaluation 

A  model  validation  study  with  respect  to  NOAA  SST 
and  NWP  from  each  NWP  product  is  presented  using 
the  statistical  metrics  explained  in  Section  4.1.  Our 
purpose  is  to  determine  whether  or  not  HYCOM  SST 


compares  with  the  NOAA  SST  climatology  better  than 
Ta  climatologies  from  each  NWP  product.  The  original 
monthly  mean  NOAA  SST  climatology  (1°)  was 
interpolated  to  the  global  HYCOM  grid  (0.72°)  for 
comparisons  with  the  model  SSTs.  Similarly,  near¬ 
surface  Ta  from  NWP  products  (ERA- 15,  ERA-40  and 
CORE-CNY)  is  also  interpolated  to  the  model  grid. 

The  mean  HYCOM  SST  error,  with  respect  to  the 
NOAA  SST  climatology  (Fig.  4a)  and  climatologies 
(Fig.  4b)  presents  striking  differences.  In  general. 


(a)  Bias:  HYCOM  SST  -  NOAA  SST 


(b)  Bias:  HYCOM  SST  -  NWP  Ta 


Fig.  4.  (a)  Annual  mean  error  (bias)  obtained  from  HYCOM  simulations  using  elimatologieal  atmospheric  forcing  from  ERA- 15,  ERA-40  and 
CORE-CNY,  when  the  model  SST  is  compared  to  (a)  NOAA  SST  climatology  and  (b)  near-surface  air  temperature  climatology  from  each  NWP 
product.  Atmospherically-forced  HYCOM  simulations  include  no  assimilation  of  any  SST  or  air  temperature.  In  the  model  simulations,  there  is  no 
relaxation  to  any  SST  or  air  temperature  climatology.  We  evaluate  time  series  of  HYCOM  SST  versus  NOAA  SST  (similarly  HYCOM  SST  versus 
NWP  T2)  climatologies  from  January  to  December  at  each  model  grid  point,  producing  the  annual  mean  bias  maps  over  the  global  ocean.  In  the  maps, 
the  white  color  represents  bias  values  between  -0.25  °C  and  0.25  °C. 
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most  of  the  mean  errors  shown  in  Fig.  4  do  not  look  like 
the  T^-T^  fields,  where  is  the  SST  used  by  the 
NWP  products,  presented  in  Fig.  1  and  discussed  in 
Section  2.2. 

To  better  visualize  whether  or  not  HYCOM  mean  SST 
error  patterns  resemble  Fa- SST  patterns  for  the  given  re¬ 
analysis  product,  we  calculate  zonal  averages  of  both 
variables  (Fig.  5).  Zonally-averaged  Fa- SST  values  can 
also  be  quite  different  for  each  of  these  re-analysis  products. 
It  should  be  emphasized  that  HYCOM  simulations  using 
ERA- 15,  ERA-40  and  CORE-CNY  forcing  have  different 
bias  values,  being  warmest  (Table  2)  with  a  globally 
averaged  bias  (HYCOM-NOAA  SST)  value  of  0.57  °C  for 
the  CORE-CNY  forced  simulation.  In  general,  the  annual 
mean  SST  bias  between  HYCOM  and  the  NOAA 
climatology  is  small  (<0.5  °C)  for  all  three  simulations 
using  the  different  atmospheric  forcing  products  over  most 
of  global  ocean. 

Most  importantly,  there  is  no  clear  relationship 
between  HYCOM  SST  bias  (i.e.,  HYCOM  SST-NOAA 
SST)  and  Fa- SST  fields  (Fig.  6)  when  using  any  of  the 
atmospheric  re-analysis  products  for  forcing  the  ocean 


model.  This  is  tme  despite  the  fact  that  the  HYCOM  SST 
bias  and  Ta-SST  values  from  a  particular  re-analysis 
product  (e.g.,  ERA- 15)  can  be  quite  different  than  those 
from  the  other  product  (e.g.,  ERA-40  and  CORE-CNY). 

If  the  model  SST  is  constrained  by  from  any 
particular  NWP  product,  then  the  differences  in  near¬ 
surface  air  temperatures  between  two  NWP  products 
should  be  similar  to  differences  in  model  SSTs,  obtained 
using  the  corresponding  NWP  forcing  in  simulations. 
This  is  certainly  not  the  case.  As  seen  from  Fig.  7,  the 
scatter  diagram  of  differences  in  between  ERA-40 
and  CORE-CNY  versus  those  in  SST  from  a  HYCOM 
simulation  forced  with  ERA-40  and  CORE-CNY  are 
quite  different.  While  differences  in  are  generally 
between  0  °C  and  1  °C,  differences  in  SST  are  generally 
between  -1.5°  and  1  °C.  There  is  only  a  weak 
correlation  (7? =-0.36)  between  the  two.  A  similar 
analysis  for  ERA- 15  versus  CORE-CNY  is  not  shown 
because  from  ERA-40  and  ERA- 15  are  almost  same 
over  the  global  ocean. 

As  further  evidence  that  bulk  parameterizations  are 
not  excessively  tracking  T^,  the  skill  of  monthly  mean 


—  HYCOM  SST  -  NOAA  SST 
_  NWP  T  -  HYCOM  SST 

a 


Latitude  belt 

Fig.  5.  Zonal  averages  of  mean  error  between  HYCOM  and  the  NOAA  SST  climatology  when  the  model  uses  climatological  atmospheric  forcing 
(wind  and  thermal  forcing)  from  the  three  NWP  products:  ERA- 15,  ERA-40  and  CORE-CNY.  Also  included  are  zonal  averages  of  SST  for  each 
product.  Zonal  averages  are  calculated  at  each  1°  latitude  belt  over  the  global  ocean. 
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Table  2 

Global  averages  of  statistical  metrics  calculated  over  the  seasonal  cycle 
between  HYCOM  and  NOAA  SST  (see  text  for  details) 


ERA- 15 

ERA-40 

CORE-CNY 

NOAA  SST 

ME  (°C) 

0.23 

0.42 

0.57 

RMS  (°C) 

0.81 

0.82 

0.98 

^cond 

0.08 

0.07 

0.07 

-^uncond 

0.29 

0.29 

0.42 

R 

0.92 

0.92 

0.92 

SS 

0.48 

0.49 

0.36 

NWP  Ta 

ME  (°C) 

1.54 

1.73 

2.32 

RMS  (°C) 

1.94 

2.08 

2.76 

B 

cond 

0.05 

0.05 

0.07 

-^uncond 

0.46 

0.54 

0.73 

R 

0.87 

0.86 

0.81 

SS 

0.25 

0.15 

-0.14 

Results  are  given  for  HYCOM  simulations  using  atmospheric  forcing 
(wind  and  thermal  forcing)  from  ERA- 15,  ERA-40  and  CORE-CNY 
Basin-averaged  means  are  calculated  over  the  entire  ice-free  global 
HYCOM  domain.  Similar  comparisons  are  made  between  HYCOM 
SST  and  NWP  air  temperature  (a  different  for  each  product).  ASS 
value  of  1  indicates  perfect  SST  match  with  respect  to  NOAA  SST  (or 
NWP  Ta)  climatology.  All  R  values  are  statistically  significant  in 
comparison  to  0.70  at  a  95%  confidence  interval. 

model  SST  against  monthly  NOAA  SST  climatology 
and  against  the  monthly  mean  re-analysis  (a  different 
Ta  for  each  simulation)  is  calculated  (Fig.  8).  We 
evaluate  time  series  of  HYCOM  SST  versus  NOAA 
SST  (similarly  HYCOM  SST  versus  NWP  air  tempera¬ 
ture)  climatologies  from  January  to  December  at  each 
model  grid  point  (see  Section  4.1).  Thus  we  produce 
spatial  variation  of  the  skill  score  over  the  global  ocean. 
In  the  panels,  blue  (red)  color  shows  negative  (positive) 
skill  scores,  and  white  color  denotes  skill  score  values 
close  to  0.  In  general,  the  success  of  HYCOM  to  predict 
monthly  SST  on  climatological  time  scales  is  evident 
when  using  atmospheric  forcing  from  any  one  of  the 
NWP  products,  in  that  there  is  positive  skill  over  most  of 
the  global  ocean.  The  SST  skill  is  even  close  to  perfect 
(i.e.,  1)  everywhere  except  the  equatorial  Pacific  and 
high  southern  latitudes.  The  low  skill  is  due  in  part  to  the 
atmospheric  forcing  and  in  part  to  deficiencies  in  the 
model,  including  relatively  coarse  model  resolution. 

The  skill  is  much  higher  against  observed  SST  than 
against  the  applied  over  the  most  of  global  ocean 
(Fig.  8b).  This  result  is  especially  evident  in  the  Indian 
Ocean  and  Atlantic  Ocean  for  the  simulations  using  the 
different  atmospheric  forcing  products.  When  HYCOM 
SST  is  validated  against  NOAA  SST  (NWP  Tf),  global 
average  of  SS  values  are  0.48  (0.25),  0.49  (0.15)  and 
0.36  (-0.14)  for  the  ERA-15,  ERA-40  and  CORE-CNY 


forcing  cases,  respectively  (Table  2).  The  reduction  in 
SS  values  is  generally  >50%  when  HYCOM  SST  is 
evaluated  with  respect  to  NPW  as  opposed  to  NOAA 
SST. 

The  low  skill  between  HYCOM  SST  and  NWP  is 
due  mainly  to  the  fact  that  ^uncond  is  significantly 
increased  (not  shown),  i.e.,  there  are  large  biases  in 
means  of  HYCOM  SST  and  NWP  based  on  the 
definition  of  unconditional  bias  (Section  4.1).  Note  that 
there  is  not  much  change  in  global  averages  of  R  and 
^cond-  Further  comparisons  of  HYCOM  SST  versus 
both  NOAA  SST  and  NWP  clearly  demonstrate  that 
the  shape  of  the  seasonal  cycle  between  two  variables 
follow  each  other  quite  well  as  evident  from  R  values 
>0.8  (Fig.  9).  The  exception  is  the  HYCOM  simulation 
with  SST  compared  to  from  CORE-CNY.  As 
expected,  all  NWP  models  have  different  boundary 
layer  parameterizations.  In  addition,  data  assimilation 
methods  and  data  type  (e.g.,  satellite  and  in-situ)  used  in 
the  assimilations  may  differ  from  one  NWP  center  to 
another  one.  Therefore,  differences  in  their  outputs,  such 
as  the  near-surface  air  temperatures  are  expected. 

Finally,  overall  HYCOM  performance  in  simulating 
climatological  SST  is  discussed  when  the  model  uses 
atmospheric  forcing  from  ERA- 15,  ERA-40  and  CORE- 
CNY.  Table  2  provides  a  summary  of  statistical  metrics 
between  HYCOM  SST  and  NOAA  SST  in  the  form  of 
global  averages.  While  the  SST  simulated  by  HYCOM 
using  atmospheric  forcing  from  CORE-CNY  did  not  yield 
results  as  good  as  those  using  ERA- 15  and  ERA-40, 
HYCOM  generally  gave  similar  success  when  using  any 
of  the  atmospheric  forcing  products. 

Although  there  have  been  many  improvements  in  the 
ERA-40  forcing  in  comparison  to  the  ERA- 15  forcing 
(Section  3.2),  the  model  response  in  simulating 
climatological  SST  remains  similar  with  global  mean 
RMS  SST  differences  of  0.81  °C  and  0.82  °C,  and  SS 
values  of  0.48  and  0.49  for  the  ERA-15  and  ERA-40 
cases.  There  is  a  warm  model  SST  bias  in  all  of  the 
simulations.  This  could  be  due  to  known  shortcomings 
in  NWP  solar  radiation  fields,  e.g.,  due  to  inaccurate 
cloudiness,  or  could  result  from  the  ocean  model 
climatology.  In  any  case,  such  a  bias  can  be  removed 
by  applying  a  spatially  varying  but  constant  in  time 
correction  to  the  total  heat  fluxes  (not  discussed  here). 

5.  Conclusions 

The  most  significant  result  of  this  paper  is  that 
virtually  all  applications  of  the  bulk  formulae,  including 
a  fixed  air  temperature,  does  not  make  the  sensible  and 
latent  heat  behave  as  though  the  model  SST  tracks  the 
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HYCOM  SST  -  NOAA  SST  (‘'C) 
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Fig.  6.  Scatter  plots  for  HYCOM  SST  bias  (HYCOM  SST-NOAA  SST)  versus  Ta-  SST  based  on  zonally-averaged  values  shown  in  Fig.  5.  The  least 
square  line  for  each  product  is  also  shown  along  with  the  linear  R  value  between  the  two  variables  for  each  product,  respectively.  The  least  squares 
lines  are  based  on  zonally-averaged  values  for  each  NWP  product.  Zonally-averaged  values  are  used  for  plotting  purposes  to  reduce  number  of  data 
over  the  global  ocean.  The  regions  where  ice  is  located  are  not  used  in  the  analyses. 
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Fig.  7.  Scatter  plots  for  SST  differences  from  two  HYCOM  simulations  forced  with  ERA-40  and  CORE-CNY  versus  differences  in  between  ERA- 
40  and  CORE-CNY.  Results  are  based  on  annual  mean  values  over  the  global  ocean  except  iee-covered  regions.  The  least  squares  line  having  a  linear 
R  value  of -0.36  is  also  shown. 
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(a)  Skill:  HYCOM  vs  NOAA  SST 


(b)  Skill:  HYCOM  SST  vs  NWP  Ta 


Fig.  8.  Dimensionless  skill  score  calculated  over  the  seasonal  cycle  for  HYCOM  SST  with  respect  to  (a)  NOAA  SST  climatology  and  (b)  near-surface 
air  temperature  climatology  from  each  NWP  product.  Results  are  shown  when  the  model  is  forced  with  NWP  products  (ERA-1 5,  ERA-40  and  CORE- 
CNY),  separately.  Note  that  air  temperature  climatology  used  for  comparisons  in  (b)  is  different  for  each  panel.  As  in  Fig.  4,  atmospherically-forced 
HYCOM  simulations  include  no  assimilation  of  any  SST  or  air  temperature,  also  no  relaxation  to  any  SST  or  air  temperature  climatology. 


near-surface  air  temperature  too  strongly.  This  is 
demonstrated  using  an  atmospherically-forced  OGCM 
(HYCOM)  with  no  assimilation  of  any  SST  and  no 
relaxation  to  any  SST  climatology. 

Accurate  SST  simulations  in  an  OGCM  depend  on 
how  the  local  changes  in  SST  made  by  the  bulk  formula 
are  modified  by  vertical  mixing  (mixed  layer  physics) 
and  advection.  Atmospherically-forced  HYCOM  simu¬ 
lations  clearly  reveal  that  while  the  RMS  SST  difference 
is  quite  small  between  SST  and  near-surface  air 
temperature  over  most  of  the  global  ocean,  the  use  of 


the  latter  in  a  bulk  parameterization  does  not  introduce 
an  inappropriate  restoring  toward  observed  SST  in  the 
model  simulations. 

Through  extensive  model-data  comparisons  using 
various  statistical  metrics,  we  have  demonstrated  that 
bulk  parameterizations  track  observed  SST  rather  than 
near-surface  air  temperature  in  atmospherically-forced 
OGCM  simulations.  When  monthly  mean  SST  simul¬ 
ated  by  HYCOM  is  compared  against  the  monthly 
mean  observed  SST  climatology  from  NOAA  and  Ta 
climatologies  from  ERA- 15,  ERA-40  and  CORE-CNY, 
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(a)  Correlation:  HYCOM  vs  NOAA  SST 
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(b)  Correlation:  HYCOM  SST  vs  NWP  Ta 


Fig.  9.  The  same  as  Fig.  8  but  for  the  correlation  coefficient.  In  the  maps,  the  white  color  represents  correlation  values  of  >0.95. 


the  model  skill  is  higher  against  observed  SST  than 
against  the  applied  used  in  the  simulations.  Overall,  the 
global  mean  of  RMS  SST  difference  for  HYCOM  is 
0.81  °C,  0.82  °C  and  0.98  °C  with  respect  to  the  NOAA 
climatology  over  the  seasonal  cycle  for  the  simulations 
using  atmospheric  wind  and  thermal  forcing  from  ERA- 
15,  ERA-40  and  CORE-CNY,  respectively.  These  global 
mean  RMS  differences  increase  by  ~150%  (1.94  °C, 
2.08  °C  and  2.76  °C,  respectively)  when  HYCOM  SST  is 
evaluated  with  against  near-surface  air  temperature  from 
the  three  NWP  products.  Model  SST  validation  against 
NOAA  SST  and  NWP  yields  similar  correlation 
coefficients  (generally  >0.90)  over  the  most  of  global 
ocean,  which  may  imply  an  indirect  relaxation  of  model 
SST  to  NWP  Ta.  However,  significantly  large  differences 


in  dimensionless  skill  score  and  RMS  differences 
between  the  pairs  of  HYCOM  SST  versus  NOAA  SST 
and  model  SST  versus  NWP  indicate  that  such  a  direct 
relaxation  is  not  the  case.  Further,  the  correlation 
between  7;(ERA-40)-7;(CORE-CNY)  and  HYCOM 
SST(ERA-40)-HYCOM  SST(CORE-CNY)  is  R  = 
-0.36  when  SST  relaxation  to  NWP  near-surface  air 
temperature  would  imply  a  substantial  positive  value. 
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Appendix  A.  New  features  of  HYCOM 

HYCOM  development  has  continued  since  the  first 
release  discussed  in  Bleck  (2002),  and  new  features 
have  been  added  to  the  model.  A  short  summary  of  the 
new  improvements  in  the  model  code  (version  2. 1 . 1 8)  is 
provided  here. 

A  hybrid  vertical  coordinate  grid  generator  (called 
hybgen)  is  used  in  choosing  the  optimal  vertical 
coordinate  at  each  location  every  time  step.  It  then 
remaps  the  vertical  coordinate  accordingly.  There  are 
significant  improvements  to  hybgen  in  the  latest  version 
of  the  model.  The  remapper  in  the  original  model 
assumed  that  each  field  was  constant  in  the  vertical 
within  each  layer.  When  remapping  layers  that  are  far 
from  the  target  isopycnal,  this  approach  can  lead  to 
excessive  diffusion.  The  current  modified  remapper  in 
HYCOM  allows  the  profile  to  vary  linearly  across  a 
layer  if  the  layer  is  not  close  to  being  isopycnal,  which 
significantly  reduces  diffusion.  In  such  cases,  the 
original  remapper  used  donor-cell  upwind  advection, 
and  the  latest  remapper  uses  the  piecewise  linear  method 
(van  Leer  1977;  Lin  et  al.,  1994)  with  a  monotonized 
central-difference  limiter.  The  hybrid  grid  generator  and 
horizontal  advection  terms  can  each  conserve  potential 
temperature  and  salinity,  potential  density  and  salinity  or 
potential  density  and  potential  temperature.  Because  the 
hybrid  grid  generator  must  act  in  density  space,  potential 
density  and  salinity  are  typically  conserved  in  both  the 
grid  generator  and  advection. 

Since  the  first  release  of  HYCOM  (Bleck,  2002), 
there  is  improved  support  for  z  and  cr  levels  in  shallow 
water,  and  more  flexible  selection  of  Laplacian  and 
biharmonic  diffusion.  There  are  several  scalar  advection 
techniques  available  in  HYCOM,  such  as  the  donor  cell, 
the  Flux-Corrected  Transport  (FCT)  scheme  (2nd  and 
4th  order)  and  Multidimensional  Positive  Definite 
Advection  Transport  Algorithm  (MPDATA)  (e.g., 
Smolarkiewicz  and  Margolin,  1998;  Balsara  and  Shu, 
2000).  While  we  do  not  discuss  it  here,  HYCOM 
simulations  demonstrate  that  FCT  is  more  computation¬ 


ally  efficient  and  less  diffusive  than  MPDATA, 
consistent  with  results  by  Smolarkiewicz  (1984).  Multi¬ 
ple  tracers  and  off-line  one-way  nesting  are  possible  in 
the  model.  HYCOM  includes  a  hybrid  to  fixed  vertical 
grid  remapper,  allowing  fixed-coordinate  nests  inside 
hybrid  coordinate  outer  domains.  Vertical  remapping 
uses  the  Piecewise  Linear  Method  (PLM)  for  fixed- 
coordinate  layers.  Isopycnal  layer  target  densities  can 
vary  spatially,  and  stability  of  these  layers  is  controlled 
by  locally-referenced  potential  density. 

HYCOM  can  be  run  using  any  one  of  five  mixed 
layer  models:  (1)  The  K-Profile  Parameterization  (KPP) 
level  1  turbulence  closure  (Large  et  al.,  1997),  (2) 
Kraus-Tumer  (KT)  mixed  layer  model  (Kraus  and 
Turner,  1967),  (3)  Price,  Weller  and  Pinkel  (PWP) 
mixed  layer  model  (Price  et  al.,  1986),  (4)  the  Mellor- 
Yamada  (MY)  level  2.5  turbulence  closure  (Mellor  and 
Yamada,  1982),  and  (5)  the  NASA  Goddard  Institute  for 
Space  Studies  (GISS)  level  2  turbulence  closure  (Canuto 
et  al.,  2002).  The  model  can  also  be  run  with  no  mixed 
layer  model. 

Earlier  HYCOM  simulations  used  the  bulk  parameter- 
izations  presented  in  Kara  et  al.  (2002),  which  has  a  few 
shortcomings  in  the  determination  of  exchange  coeffi¬ 
cients  of  sensible  and  latent  heat  fluxes  (Cg  and  Ci)  at  high 
and  low  wind  speeds.  This  was  due  to  limitations  in  the 
observation-based  input  data  sets  used  for  deriving  the 
exchange  coefficients.  Also,  the  effect  of  humidity  was 
not  included  in  the  parameterizations  of  air-sea  stability. 
With  the  availability  of  the  improved  COARE  algorithm 
version  3.0  (Fairall  et  al.,  2003),  Kara  et  al.  (2005d) 
derived  new  polynomials  for  C\  and  Cg  for  use  in 
HYCOM,  including  the  full  effect  of  stability  for  a  wide 
range  of  conditions  occurring  over  the  global  ocean.  The 
new  data  set  from  COARE  3.0  reduces  the  under¬ 
estimation  or  overestimation  in  the  polynomials  (~  10  to 
20%)  used  for  deriving  the  exchange  coefficients 
presented  in  Kara  et  al.  (2000,  2002).  In  addition,  they 
have  the  advantage  of  providing  reliable  transfer 
coefficients  at  low  and  high  wind  speeds.  They  represent 
only  an  approximation  to  the  COARE  algorithm,  but  with 
the  advantage  of  robustness  and  computational  efficiency 
that  make  it  suitable  for  use  in  various  air-sea  interaction 
applications  and  in  any  OGCM. 
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